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We experimentally study the critical properties of the non-equilibrium solid-liquid-like transition 
that takes place in vibrated granular matter. The critical dynamics is characterized by the coupling 
of the density field with the bond-orientational order parameter Q 4 , which measures the degree of 
local crystallization. Two setups are compared, which present the transition at different critical 
accelerations as a a result of modifying the energy dissipation parameters. In both setups five 
independent critical exponents are measured, associated to different properties of Qn : the correlation 
length, relaxation time, vanishing wavenumber limit (static susceptibility), the hydrodynamic regime 
of the pair correlation function, and the amplitude of the order parameter. The respective critical 
exponents agree in both setups and are given by u± = 1, i'll = 2, 7 = 1, 77 « 0.6 — 0.67, and /3 = 1/2, 
whereas the dynamical critical exponent is z = v\\/v± — 2. The agreement on five exponents is an 
exigent test for the universality of the transition. Thus, while dissipation is strictly necessary to form 
the crystal, the path the system undergoes towards the phase separation is part of a well defined 
universality class. In fact, the local order shows critical properties while density does not. Being the 
later conserved, the appropriate model that couples both is model C in the Hohenberg and Halperin 
classification. The measured exponents are in accord with the non-equilibrium extension to model 
C if we assume that a, the exponent associated in equilibrium to the specific heat divergence but 
with no counterpart in this non-equilibrium experiment, vanishes. 

PACS numbers: 64.60.Ht, 64.70.qj 05.40.-a, 45.70.-n, 


I. INTRODUCTION 

Thermodynamics and equilibrium statistical mechan¬ 
ics have shown an extraordinary success in describing 
phase transitions. It has been possible to classify them, 
predict their properties, and understand the critical phe¬ 
nomena. Also, they provide a framework that can be 
used to build efficient numerical tools as the Monte Carlo 
simulations or to design and analyze experiments, for 
example exploiting the concept of universality classes. 
From a theoretical viewpoint the critical phase transi¬ 
tions were classified in universality classes in terms of 
symmetry and conservation properties mm- 

In several order-disorder phase transitions the dynam¬ 
ics of the transition needs the interplay of two or more 
order parameters mg. One particularly interesting 
case is the so called model C in the Hohenberg and 
Halperin classification, where a non-conserved critical or¬ 
der parameter couples to the conserved noncritical den¬ 
sity [n [g m- Examples where this model can be ap¬ 
plied are varied, comprising liquid-liquid phase transi¬ 
tions Enzusmu, binary alloys [6], and anisotropic mag¬ 
nets fl2| . 

In out of equilibrium conditions there is no such sys¬ 
tematic framework that can be used to analyze or classify 
phase transitions. There are, nevertheless, some notori¬ 
ous examples of prototype models that have been shown 
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to be quite general, allowing other systems to be com¬ 
pared with them. Some of these are, for example, the 
directed percolation process, the Kardar-Parisi-Zhang 
model for surface growth and the Swift-Hohenberg model 
mm- The critical phenomena methodology in equilib¬ 
rium has been extended to these cases, where again we 
find critical exponents, universality classes, and scaling 
functions, with the dynamic renormalization group as a 
useful approach. But, still much understanding is needed 
to reach a state of knowledge comparable to the equilib¬ 
rium case. 

Granular matter, due to the strong dissipation present 
at the particle interactions and the corresponding need 
for continuous energy injection to sustain dynamic states, 
is an excellent candidate for studying out of equilibrium 
phase transitions. Important progress has been made 
along this direction [HHHZ], but a wider view in the con¬ 
text of dynamical phase transitions is still lacking. Re¬ 
cently, we presented an experimental study of a granular 
liquid-solid-like phase transition in a vibrated quasi-two- 
dimensional granular system [28] . There, we showed that 
the transition is characterized not only by the density 
field but also by a bond-orientational order parameter, 
which is described by the model C. To our understand¬ 
ing, this is the first non-equilibrium realization of this 
class of phase transition. 

In this manuscript we extend the results found previ¬ 
ously by focusing now on the universality of the transi¬ 
tion. That is, we aim to verify if by varying some ex¬ 
perimental parameters, the same critical exponents are 
found. We remark that in Ref. [ 25 ] five independent 
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critical exponents were found which, therefore, put an 
exigent test the universality condition. 

The article is organized as follows. Section [TT| describes 
the liquid-solid-like transition that takes place in con¬ 
fined quasi two dimensional granular systems. In Sec. 
|III| the experimental setup and procedures are explained, 
describing the two configurations that are used to test the 
universality of the critical exponents. The experimental 
results and the determination of the critical accelerations 
and exponents are presented in Sec. m Finally, conclu¬ 
sions are given in Sec. |Vj 


II. THE LIQUID-SOLID-LIKE TRANSITION IN 
QUASI TWO DIMENSIONAL GRANULAR 
SYSTEMS 



Granular matter is ubiquitous in our daily life, still 
a fully understanding of its dynamical behavior has re¬ 
mained rather elusive for several years [29} ]30| . Dry gran¬ 
ular matter is a collection of athermal macroscopic par¬ 
ticles that interact mainly through hard core-like dissi¬ 
pative collisions, and depending on external conditions 
such as pressure or packing fraction it may behave as a 
solid, a liquid or even a gas. 

Energy injection is needed to compensate the energy 
dissipated in grain-grain and grain-wall collisions. Vibra¬ 
tions are an efficient method to perform this task in a dis¬ 
tributed and controlled way, being possible to reach sta¬ 
ble states. Among the different possibilities to perform 
the vibrations, the quasi-two dimensional (Q2D) geom¬ 
etry has gained attention because it allows the granular 
system to be followed in detail at the individual grain 
scale together with the collective dynamics. In Q2D sys¬ 
tems, grains are placed in shallow box, with a height 
that is between one and two particle diameters, while it 
is large in the horizontal directions. The box is vertically 
vibrated with maximal accelerations larger than gravity 
such that grains acquire vertical energy that is trans¬ 
ferred via collisions to the horizontal degrees of freedom. 
In the reduced two-dimensional dynamics, the system re¬ 
sembles a liquid but with the important difference that 
the collisions are not conservative. 

For fixed density and vibration frequency, by increas¬ 
ing the maximum acceleration the system presents a 
phase transition where grains cluster in solid-like regions 
with crystalline order Ding. The instability is origi¬ 
nated in the development of an effective negative com¬ 
pressibility [261131] and the transient dynamics is gov¬ 
erned by waves [26]. Once the system is segregated 
the solid-like cluster show fluctuations in the interface 
that are well described by the capillary theory, allowing 
us to extract an effective surface tension that has non- 
equilibrium origin 152] , 

In our previous work we studied experimentally the 
solid-liquid-like phase transition in the vibrated Q2D ge¬ 
ometry. Depending on the filling fraction and height of 
the cell, the transition was either discontinuous or contin¬ 


FIG. 1. Schematic of the experimental setup, (a) Top view of 
the quasi-2D cell, with L x = L y = lOOd. (b) Side view of the 
setup. The vertical height in the cell is L z = 1.94d ± 0.02d. 
The cell is illuminated from below with a 2D array of light 
emitting diodes, which light is diffused with a white acrylic 
sheet placed between the array and the cell. (1) camera, (2) 
quasi-2D cell, (3) electromechanical shaker, (4) accelerometer. 


uous. In the later case critical phenomena develop, with 
static and dynamical properties that show power law be¬ 
havior in terms of the distance to the critical acceleration 



III. EXPERIMENTAL SETUP AND 
PROCEDURES 

In this paper we extend the analysis of the previously 
published study [25] . Two sets of experiments are used 
to test the universality hypothesis. They differ in the 
dissipation coefficients, which will be labeled experiments 
A and B, with larger and lower dissipation respectively. 

The granular system is composed of N ~ 10 4 stainless 
steel spherical particles of mass m = 4.45 x 1CU 3 g, and 
d = 1 mm diameter. The Q2D box has lateral dimensions 
L x = L y = L = lOOd. The box consists of two 10 mm 
thick glass plates separated by a square metallic frame. 
Each inner glass surface has an indium tin oxide (ITO) 
coating, which dissipates electrostatic charges generated 
by collisions of particles with the walls. The box is fixed 
to a base by four posts placed at each corner of the cell. 
The base supports an array of high intensity light emit¬ 
ting diodes. A piezoelectric accelerometer is fixed to the 
base, allowing the measurement of the imposed forcing 
acceleration with a resolution of O.Olg. The main advan¬ 
tage of this setup is that particles are illuminated from 
below. They are then visualized as dark disks on a white 
background. This allows to detect particles in dense clus¬ 
ters, even when particles are partially mounted on top of 
each other. A scheme of the setup is shown in Fig. [T] A 
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FIG. 2. A typical raw image of the complete system, for type 
B experiment and T = 4.05 ±0.01. Particles are visualized as 
black disks on a white background. 


typical image of the system is shown in Fig. [5] 

The whole setup is forced sinusoidally with an elec¬ 
tromechanical shaker, with vertical displacement z(t) = 
Asin(wt). Top view images are obtained with a camera 
at 10 fps. The images acquired have a typical resolu¬ 
tion of 1020 x 1020 pix 2 , thus we obtain particles of 10 
pix diameter approximately. Particle positions are de¬ 
termined at sub-pixel accuracy. Results have been ob¬ 
tained by fixing the particle number N, cell height L z 
and driving frequency / = ui/2n = 1/T = 80 Hz. The 
dimensionless acceleration T = Auj 2 /g is varied in the 
range 1 — 6. In this present work we focus on the config¬ 
uration C2 of reference [25], namely L z = 1.94d ± 0.02d 
and N = 11504, which implies the filling fraction to be 
<j> = Nird 2 /4L 2 = 0.904. This filling fraction value is 
possible because the system is Q2D, which allows to ac¬ 
commodate more particles than for a pure 2D system. 

As was demonstrated in [28], under this particular 
configuration the system undergoes a second-order type 
phase transition. The aim of this paper is to present a 
more detailed study of the transition itself. Addition¬ 
ally, we are particularly interested if there is universality 
in this granular out-of-equilibrium phase transition. In 
order to study the universal behavior of the critical expo¬ 
nents, we change the dissipation parameters of the system 
by changing both the bottom and top lids. We use two 
pairs of ITO (indium tin oxide) coated glass plates. In the 
experiments of type A the coating is 25 nrn thick, whereas 
for type B experiments it is 750 nrn thick. It has been 
shown that ITO films mechanical, fracture, morphologi¬ 


cal and electric properties depend on fabrication details 
such as thickness, annealing, substrate, feed gas, among 
others [38]. So, it is expected that the dissipation param¬ 
eters, namely inelastic and friction coefficients, should be 
different for the two sets of ITO coated glass plates. As 
will be shown below, our results indeed demonstrate that 
the dissipation is stronger for the thinner ITO coated 
glass plates (type A experiments). 

The ITO coating works very well for many hours of 
experimental runs. Eventually, it does however get dam¬ 
aged by particle collisions. The precise time of repro¬ 
ducible runs is probably function of the ITO coating 
thickness. This is supported by the thickness dependence 
of the crack onset strain in ITO thin films, which is lower 
for thicker coatings [33] ■ All data presented in this paper 
correspond to reproducible runs for which no important 
damage was noticeable. In fact, a surface damaged ITO 
coating is manifested in important changes of the mea¬ 
sured quantities - such as the density and order structure 
factors - with respect to those obtained for a new pair of 
ITO coated glass plates. We conjecture that the damage 
occurs because of erosion of the ITO coating, which in 
turn affects particle interactions by an increase of elec¬ 
trostatic forces and contamination of the system by dust 
formed from the ITO coating. In order to insure repro¬ 
ducibility, glass plates were changed periodically during 
the time duration of all the experimental runs, and all 
parts of the experiments are cleaned in an ultrasonic bath 
before mounting the experimental cell again, including 
the particles. 

Two important issues are the top and bottom plates 
flatness and homogeneity of L z through the experimental 
cell. With the cell empty of particles the height is mea¬ 
sured at nine different positions on a regular spaced grid 
with an optical microscope (see supplementary material 
of [28]). Within experimental errors the vertical height 
is homogeneous, L z = 1.94d ± 0.02d. However, the ho¬ 
mogeneity of L z does not insure that the glass plates 
are both flat. Because of the mechanical constrains and 
stresses exerted on the plates, some small residual curva¬ 
ture could exist. 

In Fig. [3](a) we show an average of 3270 images for 
r = 2.05, well below the solid-liquid transition (for type 
B experiments, F c ss 4.5). Only a small part of the setup 
is shown, of size 15 d x 15d. Viewed from above it corre¬ 
sponds to the left side wall at almost half of the cell. It 
is clear that near the side wall particles tend to form lay¬ 
ers of thickness ~ d, within which particles barely move 
compared with the rest of the system. This layering, 
observed at all side walls, is induced by the extra dissi¬ 
pation given by the side wall collisions, similar to what 
has been shown in sheared granular matter [23]. How¬ 
ever, for a distance > 5 d from a side wall, the system 
is nearly homogeneous. The inset of Fig. |3][a) shows an 
average intensity plot as a function of the position (time 
and y averaged), showing the layers close to the edge and 
then saturating at a constant value. Also, by means of a 
time average coarse-graining procedure we can compute 
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FIG. 3. (a) Average image obtained from 3270 raw images. Only a section of the setup near a wall is shown (of size 15d x 15d). 
Regions in dark grey (low intensity) results from more particles in these zones. The vertical stripes at the left account for 
layers of particles that move less and spend more time near the wall because of the extra dissipation at the side wall. The 
inset shows the intensity of this image as a function of position, showing that for x > 5d the system is almost homogeneous, 
(b) Time averaged coarse-grained particle density obtained directly from particle detection as a function of x, with d/4 as bin 
width. It also shows oscillations due to side wall layering and confirms the homogeneity for x > 5d at this scale, (c) The same 
coarse-grained density at a different scale, showing that the density at the center is about 2% higher than close to the edges. 
Here the bin width is d/2. For both (b) and (c) Ni is the average number of particles detected in a vertical bin spanning 
y = [0, lOOd] and Xi = [(* — l)wb,iwb), where Wb = d/4 or d/2 is the bin width. No is defined as the number of particles that 
would fit into a bin with an homogeneous density p 0 = N/L 2 . 


the density of particles as a function of distance from a 
wall. This is shown in Figs. §b) and [3jc). Close to the 
wall we notice the same particle layering and an apparent 
homogenization for x > 5d (Fig. [3](b)). The asymptotic 
limit is slightly less than expected to compensate for the 
particle excess at the side walls. At a larger scale, as 
shown in Fig. [3](c), we observe that there is a small 
density gradient that leads to an excess of particles at 
the cell’s center compared to the edges: density is about 
2% higher at the center that near the side walls. This 
can be the result of three effects, which are probably all 
present. Firstly, even though the cell’s height is homoge¬ 
neous within experimental errors, there might be a small 
residual concavity, which makes the particles to accumu¬ 
late near the center. A second factor could be that the 
vertical acceleration may not be constant throughout the 
cell. In other words, there could be some “flapping” of 
the system. Indeed, we have measured the acceleration of 
the cell at different positions and we find that it is about 
0.2% higher near the borders that at the center. Finally, 
the presence of dissipative side walls could induce such 
large scale inhomogeneity. 

Another important issue is the mechanical leveling of 
the whole setup (for details, see supplementary material 
of Ref. |2B]). The cell’s horizontality, and thus the sys¬ 
tem’s isotropy, can be checked by two ways: First, it 
is verified through the static structure factor S(k ) (de¬ 
fined below) and the average density Fourier components 
( p(k,t )) (for now we can say that S(k) is a measure of 
density fluctuations in Fourier space). When these quan¬ 
tities are plotted versus k x and k y there is no preferential 


direction. For example, S(k x ,k y ) shows the characteris¬ 
tic symmetric circular annular shape at kd = 2i r, where 
k = \k\. A second verification is given by the even sym¬ 
metry of the coarse-grained density with respect to the 
cell center line, as can be observed from the data pre¬ 
sented in Fig. [3](c), which is indeed symmetric with re¬ 
spect to x/d = 50. 


The particle detection is done by using a modified 
open source Matlab code which uses a least-square al¬ 
gorithm [HI] . Our modified version in C++ & CUDA 
allows faster computation for large number of particles 
[151 35^. The algorithm allows to detect both layers of 
particles in a dense solid cluster, where the top layer par¬ 
ticles are placed in the valleys that the bottom particles 
form. Typical experimental runs consist of at least 30 
video acquisitions, one for each A, of about 3300 images 
each. Therefore, the complete number of images to ana¬ 
lyze for a single experimental run is about 10 5 . 


Finally, for fixed geometry, particle density and vibra¬ 
tion frequency we perform vibration amplitude ramps, 
from r > 1 in the liquid phase, increasing A going 
through the solid-liquid transition that is reached at 
r = T c . In order to check that the transition is continu¬ 
ous, for some runs we also perform decreasing amplitude 
ramps starting above r c . 
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FIG. 4. (Color online) Color plots of log 10 ((p(k,t)p(k,t)*)/N) (a), log 10 ({p{k, t))(p(k, t)*)/N) (b) and the structure factor 
function S(k) (c) as functions of k x and k y . From S(k) the isotropy of the fluctuations is clear from the circular halo shape of 
the pre-peak, (d) S(k) for a wide range of kd, obtained from S(k x , k y ) by simply plotting as a function of k = yjk\ + k y . For 
all these figures T = 3.75 < F c and data was obtained from type B experiment. 


IV. EXPERIMENTAL RESULTS 

A. Static structure function. 


Particle positions fj(t) in the plane (x,y) are deter¬ 
mined for each time t. Experimentally, there is no access 
to the z coordinate. Thus, the 2D microscopic density 
field Fourier components are 

r N . 
p(k,t) = / d 2 re if k p(f,t) = ( 1 ) 

3 =1 


The static structure factor S{k) measures the intensity 
of density fluctuations in Fourier space: 


S{k) 

S(k) 


(\p(k,t) - (p(k,t ))| 2 ) 

N 

(p{k,t)p(k,t)*) - (p{k,t))(p(k,t)*) 
N 


( 2 ) 

(3) 


where ( ) denotes time averaging. In general 
(p(k, t)) 7 ^ 0 , due to inhomogeneities induced by bound¬ 
ary conditions, as those shown in Fig. [3] The wave 
vectors are computed from k = n(n x i + n y ])/L , where 
n x ,n y £ N, and k = |fc|. 
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FIG. 5. (Color online) (a) S(k) in the large wavelength limit for different accelerations F. The pre-peak grows and shifts 
towards lower k for increasing F. We define the pre-peak’s maximum S ma x that is obtained at k = k*. The width A of the 
pre-peak is defined as the width at half height. (b,c) and £/d = n/k* as functions of T for both types of experiments, (d) 

A d versus d/£. A linear dependence is observed for £ > 15d (d/C 0.07). 


In Figs. |d|a) and[l](b) we present color plots of the two 
terms that are used for the computation of S(k x , k y ). 
namely ( p(k,t)p(k,t)*)/N and (p{k, t))(p(k, t)*)/N as 
functions of k x and k y (both in log 10 scale). It turns 
out that both quantities are strongly non-monotonic and 
are different by several orders of magnitude. In the lower 
wavenumber range that is plotted these quantities show 
a set of peaks placed on a regular grid on top of a smooth 
background, taking large values when both n x and n y are 
odd, whereas the other modes are much lower. This is 
understood by the even symmetry that the density field 
has with respect to the cell’s center, as shown in Fig. 


Sc). Thus, its Fourier decomposition yields that the not 
pure odd harmonics should vanish, having very low val¬ 
ues in practice. In Fig. [4^c) , the subtraction of these two 
terms is shown, which defines S(k x , k v ). The smoothness 
of the resulting function, with no discrimination between 
even or odd modes, indicates that the density fluctua¬ 
tions are governed by long wavelength dynamics and not 
by the static density profile. 

The structure factor presents a notorious pre-peak cen¬ 
tered at kd ~ 0.2, which corresponds to a large wave¬ 
length structure of size ~ 15d (the pre-peak term refers 
to the peak being at lower wavenumber than the one cor- 
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responding to the first coordination shell at kd = 27r, see 
below). The associated density fluctuations are indeed 
visible by simple visual inspection (see Fig. i- Fig- Sd) 
shows S(k) for a wider range of kd , up to kd ~ 40. It 
has the usual form expected for liquids with short range 
order, presenting a maximum at kd = 2it related to 
the first coordination shell as well as the pre-peak dis¬ 
cussed before. This pre-peak can be located in the range 
kd = 0.1 — 0.3 depending on the exact value of T. Simi¬ 
lar density fluctuations have been observed in amorphous 
materials mm, which have been consistently related to 
the existence of medium-range-crystalline-order. These 
density fluctuations have also been related to the pres¬ 
ence of a pre-peak in the structure factor. In fact, in the 
amorphous literature the pre-peak is known as the First 
Sharp Diffraction Peak (FSDP) because it appears at low 
k and is obtained from X-ray diffraction measurements. 

Fig. |5](a) presents S(k) for small wavenumbers at dif¬ 
ferent accelerations T below r c . From this figure it is 
clear how the pre-peak evolves as the system is driven 
towards the transition by increasing its acceleration to¬ 
wards r c : its height increases and its position shifts to 
lower k. This implies that density fluctuations become 
larger in size and stronger as F increases. We character¬ 
ize the pre-peak by its maximum value Smax = S(k*), 
and the associated characteristic length scale £ = tt/ k*, 
where k* is the position of the pre-peak. These quanti¬ 
ties are plotted in Figs. §b) and [5](c) as functions of T 
for increasing amplitude ramps and for both experiment 
types A and B. Although the data points are scattered 
and that they do not really overlap, specially at higher F, 
both quantities show similar trends for each experimen¬ 
tal type. Both S max and £/d increase as the transition 
is approached, although the former seems to saturate at 
larger T. We observe no great differences between the 
two ITO coatings, being their final values (near the tran¬ 
sition) very similar, S max « 0.6 — 0.1 and £/d « 20 — 35. 
We also present in Fig. [5](d) the width of the pre-peak, A, 
defined as its width at half S max . This quantity is a mea¬ 
sure of the dispersion around the characteristic length £. 
The collapse and scatter of the data are improved with 
respect to the other quantities, with no dependence on 
the different dissipation parameters. 

By observing visually the persistence of the solid clus¬ 
ters we conclude that for the A type experiment, the 
transition is located at T c ~ 5.1, whereas for type B 
it is found to be F c ~ 4.5. Given that we are dealing 
with a continuous phase transition these are just approx¬ 
imated values. However, neither S mlLX nor £ show evident 
changes at these values. 

As a conclusion to this first part we can say that den¬ 
sity fluctuations do not show critical behavior, but they 
are needed to create regions of high order. This is also 
evident from visual inspection; higher density patches are 
indeed more ordered, as can be verified in Fig. [2j Den¬ 
sity is a conserved field. Its fluctuations are however lim¬ 
ited by the system’s vertical geometrical constrain and 
the fact that the particles are in practice hard spheres. 
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FIG. 6. Schematic representation of two square interlaced 
layers for which particles are closed packed. The center of 
particles in the bottom (top) layer are shown with solid (open) 
black circles. The 2D projected square lattice has a unit cell 
length d/y/2 . 


In what follows, medium range order will be analyzed 
with an appropriate bond-orientational order parameter, 
which does indeed present critical behavior. 


B. Bond-orientational order parameter. 

In the vicinity of the transition, fluctuations of high 
density present the same square symmetry as the solid 
phase. In the quasi-2D geometry the solid phase con¬ 
sists of two square interlaced layers instead of the hexag¬ 
onal layer that is characteristic of 2D systems [25]. The 
local order can be characterized through a 4-fold bond- 
orientational order parameter. This is still valid in quasi- 
2D geometry because the interlaced two-layer square lat¬ 
tices (with unit cell length d in each plane) result also 
in a square lattice when projected in 2D, with unit cell 
length d/y/2 when the grains are close packed, as shown 
in Fig. [6] The 4-fold bond-orientational order parameter 
per particle is defined 

Nj 

= (4) 

3 S = 1 

where Nj is the number of nearest neighbors of particle 
j and a 3 s is the angle between the neighbor s of parti¬ 
cle j and the x axis. For a particle in a square lattice, 
IQH = 1 and the complex phase measures the square lat¬ 
tice orientation respect to the x axis. For details on the 
computation of Q\ we refer the reader to the supplemen¬ 
tary information of Ref. [25] . 

Representative maps of |Q^| are shown in Fig. [?] for 
three accelerations, T < r c , T « r c and T > r c . In 
this case the maps are obtained from images for exper¬ 
iment type A (T c ~ 5.1). Below the transition the or¬ 
dered patches, or crystallites, are first small, more or less 





FIG. 7. (Color online) Color maps of the absolute value of the 4-fold bond-orientational order parameter in real space (for each 
particle we plot IQ 4 I) for T = 4.18 (a), 5.10 (b) and 5.42 (c) (Experiment type A, r c « 5.1). In each figure only a part of the 
system is shown, from 0.2L to 0.8L in each horizontal dimension. The coloring is detailed in the legend, with the most ordered 
particles in red (\Qi\ = 0.9 — 1). Solid (open) circles correspond to particles classified as particles in the solid (liquid) phase. 
This classification is performed by measuring the area of the Voronoi area of each particle (see supplementary document of 

mr 


homogeneously distributed in space and are of short life¬ 
time. They increase in size and live for longer time as 
T approaches T c . Also, they tend appear more near the 
center than at the sidewalls, which we relate to the large 
scale small density inhomogeneity discussed before. The 
quantitative study of crystallites size and lifetime is pre¬ 
sented below. 


As discussed before, the local density can not change 
strongly at the transition because of the system’s vertical 
geometrical confinement. Thus, the correct order param¬ 
eter must be the local 4-fold symmetry order, measured 
through the orientational parameter Q\. With this in 
mind, we define its global average, in space first and then 
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FIG. 8. (Color online) (a) Global average of 4-fold orientational order parameter (IQ 4 I) versus T for the two ITO coatings: 
Experiment type A (thin ITO) with increasing T ramps (□) and experiment type B (thick ITO) with increasing (o) and 
decreasing (•) T ramps. Continuous lines show the linear trend fit for 2.5 < T < T c and the supercritical deviation for T is T c . 
The adjusted critical accelerations are T c = 5.12 ± 0.01 (type A) and T c = 4.48 ± 0.03 (type B). (b) AQ 4 = (IQ 4 I) — Q\ versus 
e = (T — r c )/r c in log-log scale for each ITO coating thickness (type A: □; type B o for both increasing and decreasing T 
ramps), where Q 4 = aT + fe is obtained from the linear trend below T c . For the thin ITO coating, a = 0.011 ± 0.001 and 
b = 0.380 ±0.002. For the thick ITO coating, a = 0.016 ±0.001 and b = 0.359 ±0.002. In sake of clarity just one representative 
error bar is shown for each case. The straight lines are power laws with exponents equal to 0.4 (dash-dotted), 1/2 (continuous) 
and 0.6 (dashed), shown as guides to the eye. 


in time, 

< |£ w> = (^f>ii), (5) 

which measures the fraction of particles in the ordered 
phase. 

In order to show that the transition is indeed of second 
order for both experiment types, we present in Fig. [ 8 jt 
the global average (IQ 4 I) versus T for the thin and thicker 
ITO coated plates (type A and B respectively). For the 
latter, both increasing and decreasing T ramps are pre¬ 
sented, showing good reproducibility. This figure indeed 
demonstrates that both configurations present a second 
order type transition, continuous and with no hysteresis. 

The qualitative behavior is the same for both experi¬ 
ment types. First, (IQ 4 I) has a linear dependence on T 
below the transition. This reflects the fact that the frac¬ 
tion of particles that form crystallites with square fold 
symmetry, even if it is transiently, increases when the 
transition is approached. For both experiments A and 
B there is a clear deviation from this linear trend above 
a given threshold. The critical acceleration -defined for 
now as the value where the qualitative change occurs- 
for the thicker ITO coating is lower (r c ~ 4.5) than the 
one for the thin ITO coating (r c ~ 5.1). Also, the initial 
linear slope below the transition is larger for the thick 
ITO coating case. Both facts, the lower critical value 


10 2 



kd 

FIG. 9. (Color online) 4-fold bond-orientational structure fac¬ 
tor S4 (/c) for several T for experiment type A (the results for 
B are basically the same). The vertical axis is in log 10 scale. 
Curves obtained for T < F c are in blue (dark gray) and for 
r > F c are in red (light gray). For both ITO thicknesses and 
for r < r c , all curves show an Ornstein-Zernike-like behavior 
in the limit kd < 1, S 4 (fc) « S 4 (0)/[1 + (&k) 2 ]. For T > F c 
the curves tend to collapse together. 
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FIG. 10. (Color online) 54 ( 0 ) (a) and £4 /d (b) versus e for experiments type A (□) and type B (o). Continuous lines are critical 
power laws, with exponents 7 = = 1 for 54 ( 0 ) and £ 4 , shown as guides to the eye. The fitted critical accelerations for 54 ( 0 ) 

and £4 /d are: F c = 5.09 ± 0.07 and F c = 5.24 ± 0.08 respectively for experiment type A; r c = 4.43 ± 0.06 and T c = 4.58 ± 0.06 
respectively for experiment type B. More details are provided in Table [I] 


and larger slope for the thicker ITO coating, are con¬ 
sistent with a lower effective dissipation at the top and 
bottom walls. Indeed, in this case the transition occurs 
at a lower amplitude, thus at lower energy injection and 
dissipation rates, and transient crystals form and grow 
more easily as T increases. In this figure, the continuous 
lines correspond to fits of a linear dependence for T < F c 
and a supercritical deviation for T > r c (details in the 
figure caption). 

The deviation from the linear trend observed for T < 
T c is defined as AQ 4 = (IQ 4 I) — Q 4 , where Q 4 is de¬ 
fined as the extrapolation of the linear trend over the 
complete range of T. Fig. §b) presents AQ 4 versus 
e = (r — r c )/F c in log-log scale for each ITO coating 
thickness. The continuous line shows the supercritical 
law AQ 4 oc Vr — T c as a guide to the eye. 

The results of Fig. [ 8 ](a) are fitted to the supercritical 
law AQ 4 = cVr — F c . The adjusted parameters are c = 
0.029T0.002 (type A) and c = 0.024T0.002 (type B). We 
conjecture that the different adjusted c values also reflect 
the difference of dissipation parameters that control the 
particle-wall collisions. 

Next, in order to characterize quantitatively the or¬ 
dered patches shown in Fig. [TJ in particular their typical 
length and time scales, we analyze the orientational order 
parameter in momentum space. Its Fourier components 
are 

N 

QS,t) = Y,Q^' Pi{t) - ( 6 ) 

1=1 

Then, local order can also be analyzed through its fluc¬ 
tuations in Fourier space by means of the 4-fold bond- 


orientational structure factor 

q ( f A (\Qi{k,t) - {Q 4 (k,t))\ 2 ) 

Si{k) = - — -. (7) 

This structure factor is shown in Fig. [9] for several ac¬ 
celerations and for experiment type A. Results for exper¬ 
iments type B are basically the same and are not shown. 
For both experiment types and for T < r c , S 4 (k) shows 
an Ornstein-Zernike-like behavior in the limit kd -C 1, 


S 4 (k) 


Si (0) 

l + (£ 4 fc) 2 ’ 


( 8 ) 


where £4 and ^ 4 ( 0 ) are the 4-fold bond-orientational cor¬ 
relation length and static susceptibility respectively. 

In Fig. 10 we present both 5 4 (0) and £ 4 /d for T < F c , 
obtained from the fits of the Ornstein-Zernike behavior 
for both ITO coatings. Both quantities are plotted as 
functions of the reduced acceleration e = (T c — r)/F e , 
where F c is obtained from a specially adapted fitting 
procedure [25]. The correlation length and susceptibil¬ 
ity vary strongly as the transition is approached. For the 
two ITO coatings these quantities obey critical power 
laws. In the limit e —> 0 they both saturate, presumably 
due to the system’s finite size. For e < 3 x 10 -2 they sat¬ 
urate to 54 ( 0 ) « 10 — 20 and £ 4 /d « 10 respectively. The 
critical-like behavior is fitted with both free and fixed 
exponents. As discussed in our previous work [25], the 
precise measurement of T c and the critical exponents is 
far from trivial. Initial fits give T c in the range 5.1 — 5.6 
and exponents that can vary up to a factor of almost 2 . 
Additionally, the fitted T c can be quite different depend¬ 
ing from which quantity they are obtained. The lack of 
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TABLE I. Adjusted parameters with free or fixed critical exponents for both ITO coatings. The fitted power laws are AQ 4 = 
c(T — T c ) /3 , S4 (0) = as r , £4 /d = bs~ I ' ± , T4/T = ce _I/ H and 84 (h) = C 00 /k 2 ~ rl . The fit of the reiaxation time 74 corresponds 
to kd = 0.09. For /3 ( 7 ) only the fixed (free) case was analyzed. 


Experiment type A - Thin ITO coating - larger wall dissipation 
Free critical exponents 

7 = 0.95 ±0.15 
T c = 4.94 ±0.22 

_ a = 0.47 ±0.03 

Fixed critical exponents 

P = 1/2 7=1 u ± = 1 i/|| = 2 

r c = 5.12 ± 0.01 r c = 5.09 ±0.07 F c = 5.24 ± 0.08 F c = 5.12 ±0.07 

c= 0.029 ±0.002 a = 0.47 ±0.01 b = 0.41 ± 0.01 5=1.10 ±0.03 

Experiment type B Thick ITO coating - lower wall dissipation 

Free critical exponents 

7 = 1.10 ±0.16 isj_ = 1.10 ±0.03 i/ii* 1.98 ±0.21 7 — 0.87 ± 0.01 (F = 4.29) 

r c = 4.54 ±0.18 F c = 4.73 ±0.03 F c = 4.55 ±0.18 7 = 0.67 ± 0.01 (F = 4.41) 

_a = 0.39 ±0.03_5 = 0.35 ±0.03_5 =1.14 ±0.07 7 = 0.63 ± 0.01 (r = 4.52) 

Fixed critical exponents 

P = 1/2 7=1 u ± = 1 i/|| = 2 

T c = 4.48 ±0.03 T c = 4.43 ±0.06 F c = 4.58 ± 0.06 F c = 4.46 ± 0.03 

c= 0.024 ±0.002 a = 0.41 ±0.01 b = 0.37 ±0.01 5 =1.17 ±0.01 


v± =0.97 ±0.21 
F c = 5.11 ±0.35 
b = 0.41 ±0.04 


i'll = 1.98 ±0.50 
T c = 5.06 ±0.38 
5= 1.11 ±0.13 


7 = 0.69 ±0.01 (F = 5.04) 
7 = 0.67 ±0.01 (F = 5.10) 
7 = 0.59 ±0.01 (F = 5.12) 


precision is due to the arbitrariness in the choice of the 
range T to be used for the fit. In [28] we present a robust 
method for the determination of these quantities, readers 
are referred to its supplementary information document 
for details. 

The results for both experiment types are given in Ta¬ 
ble |T] (the columns referring to the relaxation time 74 
and S^k) in the hydrodynamic regime will be discussed 
below). We conclude that for both ITO coatings, the 
exponents are the same within experimental errors. 

In what follows, the exponents are assumed to be fixed 
and we define 


parameters. This supports the choice of fixing the critical 
exponents. 

Now, we present the characterization of the crys¬ 
tallite’s relaxation time. Time correlations are com¬ 
puted through the two-time bond-orientational correla¬ 
tion function 


F 4 (fc,r) 


(SQ 4 (k, t ± t) 6 Q A (k, t)*) 

~lsr 


( 10 ) 


where * stands for the complex conjugate and 6 Q 4 (k, t) = 
Q 4 (k,t) — (Q 4 (fc, t)). Our results show that for low 
wavevectors 


S 4 (0) = 5e" * 7 , ii/d = b £~ v± , (9) 


F 4 (k,r) ss F 4 (k,0) exp(-r/r 4 (fc)), (11) 


with the critical exponents 7=1 and i/j_ = 1. The crit¬ 
ical divergence with e makes it necessary to fit the T 0 
separately for each case. For experiment type A the ad¬ 
justed critical accelerations are T c = 5.09 ± 0.07 and 
T c = 5.24 ± 0.08 for 5 4 (0) and £ 4 respectively, whereas 
for B they are r c = 4.43 ± 0.06 and r c = 4.58 ± 0.06 
respectively. For each ITO coating we find that within 
experimental errors both critical accelerations are very 
consistent, as well as with the value obtained from the 
supercritical behavior of A Q 4 (T c = 5.12 ± 0.01 and 
r c = 4.48 ± 0.03 for A and B respectively). Notice that 
now the critical accelerations are obtained from fits of 
measured quantities below the transition, whereas for 
AQ 4 they were obtained with fits of the order param¬ 
eter above the transition. Finally, the critical accelera¬ 
tions are less consistent between the different measured 
quantities when the critical exponents are let to be free 


from which the relaxation time r 4 (fc) is measured. The 
exponential decay behavior is shown in Fig. El From the 
measured relaxation times we also obtain a critical-like 
behavior, which is presented in Fig. [12] for the particu¬ 
lar cases kd = 0.09, 0.13 and 0.17. Here, again we obtain 
critical power law forms, for both ITO coatings. We have 
also used free and fixed critical exponents, for which de¬ 
tails are given in Table [I] The conclusion is the same as 
for the correlation length and susceptibility, namely that 
within experimental errors the free exponent adjustment 
is very consistent with a fixed critical exponent v\\ = 2 . 
Thus, we assume this critical exponent to be fixed and 
we define 

T A /T = C£- V II, (12) 

with zz|| = 2 and T is the vibration period. The adjusted 
critical accelerations are F e = 5.12±0.07 and r c = 4.46± 
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range d /£4 ^ fcd ^ 1, for configuration A, the measured 
critical exponents are 77 = 0.69±0.01, 77 = 0.67±0.01 and 
77 = 0.59 ± 0.01 for T = 5.04, 5.10 and 5.12 respectively. 
Similarly, for experiment type B they are rj = 0.87±0.01, 
T] = 0.67 ± 0.01 and 77 = 0.63 ± 0.01 for F = 4.29, 4.41 
and 4.52 respectively. Thus, although the anomalous ex¬ 
ponent ij varies rather strongly depending on T c — T, we 
can state that close enough to T c it can be bounded in 
the range 0.6 — 0.67 for both ITO coatings. In conclu¬ 
sion, there is clearly a hydrodynamic regime for which 
the power behavior is valid, even for a wider range than 
predicted. However, the measurement of 77 needs to be 
done extremely close to T c . With the present data, we 
can state that 77 « 0.6 — 0.67 is a good estimation. 


V. DISCUSSION AND CONCLUSIONS 


FIG. 11. (Color online) Dynamic 4-fold bond-orientational 
structure factor F<±(k, t ) for several kd and F = 4.36 < F c (ex¬ 
periment type B), which shows an exponential decay. Large 
wavelengths (small wavenumbers) decay slower than short 
wavelengths (large wavenumbers). The behavior of F 4 .(k,r) 
is basically the same for both ITO coatings. Only data above 
de noise level is shown for each kd. 


0.03, for A and B respectively. The relaxation time also 
seems to saturate for small e, which occurs at smaller e 
for lower fc, that is for fluctuations of larger size. 

From the relaxation time analysis we conclude that the 
dynamical exponent is z = i'\\/i , ± = 2 . As usual, there is 
critical slowing down in the dynamics. As a consequence, 
close to the critical point, stationary states are obtained 
after a long relaxation has taken place. Taken that into 
account, all T ramps are chosen to be slow. Also, averages 
are taken for long times. 

As a final evidence of the observed criticality we now 
turn to the characterization of the anomalous exponent 77 . 
In the hydrodynamic regime, d /£4 -C kd -C 1, Si(k) is 
expected to present a power law decay 

(is) 


where rj is the critical exponent related to the decay of 
the pair correlation function g(r) ~ r D ~ 2 -with D the 
dimensionality. 

Fig. 13 presents Si(k) in log-log scale for various F. 
Indeed, as the transition is approached, curves tend to 
collapse for shorter wavelengths. They are clearly dif¬ 
ferent for larger wavelengths as they converge to differ¬ 
ent static susceptibilities ^(O). In principle, 77 must be 
obtained in the limit T —> F c . As the critical acceler¬ 
ation is not known with sufficient precision, in Fig. |T3] 
we present 64 (fc) for three accelerations in the vecinity 
of r c for both ITO coatings. For the highest T, it is 
not even certain that the system is below the transition 
or not. Performing power law fits constraining fc to the 


We studied the liquid-solid-like transition that takes 
place in confined quasi two dimensional granular sys¬ 
tems. Two configurations that differ in the dissipation 
were analyzed, presenting the transition at different crit¬ 
ical accelerations. We have demonstrated that the non¬ 
equilibrium transition is a second order type for both 
configurations studied. 

Besides, we showed that in our experiments, for both 
cases, density fluctuations do not present strong varia¬ 
tions at the transition. In fact, the static structure fac¬ 
tor S(k) actually presents a peak at low wavenumbers, 
which is related to the existence of medium range crys¬ 
talline order m ■ In our case, the characteristic length 
£ of the these structures in the system does not show 
critical behavior. 

On the contrary, local order presents critical behavior. 
It is characterized through the bond-orientational order 
parameter Q4, which in Fourier space shows an Ornstein- 
Zernike-like behavior. The associated correlation length 
£4, the relaxation time T4, the zero fc limit of Q4 fluctua¬ 
tions (static susceptibility), the pair correlation function 
of Q 4, and the amplitude of the order parameter obey 
critical power laws, with saturations due to finite size 
effects. Their respective critical exponents are v j_ = 1, 
i/|| = 2, 7 = 1, 77 ~ 0.6 — 0.67 and ft = 1/2, whereas the 
dynamical critical exponent 2 = v\\/v±_ = 2. Although 
the critical accelerations and the pre-factors of the power 
laws differ between the two setups, the reported criti¬ 
cal exponents are the same. Hence, while dissipation is 
strictly necessary to form the crystal, the path the sys¬ 
tem undergoes towards the phase separation is part of a 
well defined universality class. 

In equilibrium, the scaling hypothesis predicts rela¬ 
tions among the critical exponents. It is worth mention¬ 
ing that the relation 7 = (2 — ij)u± is not satisfied, while 
a + 2(3 + 7 = 2 and u±D = 2 — a (D = 2 is the spatial 
dimension) can be satisfied simultaneously if a = 0. This 
exponent, associated in equilibrium to the specific heat 
divergence, has no interpretation out of equilibrium. 

The critical order parameter in the present case is a 
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FIG. 12. (Color online) ta/T versus e for experiment type A (left) and B (right) for three wavenumbers. The continuous 
lines are critical power laws with exponents equal to —2, shown as guides to the eye. The fitted critical accelerations are 
T c = 5.12 ± 0.07 and T c = 4.46 ± 0.03 for A and B respectively. 



FIG. 13. (Color online) ^(fc) in log-log scale for several Y for experiments type A (left) and type B (right). The continuous line 
corresponds to r\ = 0.63. We recall that F c as 5.1 and F c as 4.5 respectively. As expected, the hydrodynamic regime becomes 
wider as the transition is approached and Si(k) obeys a power law. 


non-conserved complex scalar field. Its dynamics, how¬ 
ever, is not expected to be autonomous even close to the 
critical point as density fluctuations are needed to create 
the ordered phase. Although it has been shown that the 
transition dynamics is mediated by waves (26] . momen¬ 
tum density decays fast due to friction. Therefore, the 
most appropriate description in the theory of dynamical 
critical phenomena is model C, in which a non-conserved 
order parameter is coupled to a conserved non-critical 


density [Jj. In this case in US, and in extensions to 
non-equilibrium dynamics |39| . the dynamical exponent 
is predicted to be z = 2 + a/v±, consistent with the 
measurements if a = 0. 
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